%%%%%%%%%%%%
% Generate historical moment conditions based on daily vol and return data
% HV --- historical return variance
% HO --- historical implied vol change variance
% HG --- covariance
% 3 horizons: 12,3,1 months
%
% Liuren Wu, liurenwu@gmail.com
%%%%%%%%%%%%
format compact; clearvars;

load('./data/Sampledailydata.mat');%,'IV','LFS','spot','ddate','mat','d1','pdelta','currency','nc','nm','nk','cdelta','-mat');
% a short sample of data on implied vols and currencies, 
% for illustration of the code only, do not distribute

% IV:  implied vol quotes, Txnmxnk*nc 
% T=sample length; nm=5 (maturity); nk=5 (delta); nc=2 (currencies);
% LFS: log forward spread ln(F/S), T x nm x nc, 
% spot:spot exchange rate T x nc
% ddate: date 
% mat: maturity (5) in years
% pdelta: put forward delta, decimals (10,25,50,75,90), cdelta=1-pdelta
% d1:     standardized variable d1 in Black-Scholes formula, translated from delta

T=length(ddate)
datestr(ddate([1,end]))


%Expand d1 and maturity index to match daily data
d1m=NaN(T,nm,nk,nc); matm=d1m;
for m=1:nm
    d1m(:,m,:,1)=repmat(d1(:)',T,1);d1m(:,m,:,2)=repmat(d1(:)',T,1);
end
for k=1:nk
    matm(:,:,k,1)=repmat(mat(:)',T,1);matm(:,:,k,2)=repmat(mat(:)',T,1);
end

m1d=mat-1/252; 
dls=spot*NaN;dls(2:end,:)=log(spot(2:end,:)./spot(1:end-1,:)); %log sport change over previous day


IVk=IV*NaN; %the IV of previous strike
dIVk=IV*NaN; %percentage implied vol change of previous day contract
dFm=IV*NaN; %percentage forward return of previous day contract

hh=[252,63,21]';nh=3; %moment condition horizon in business days
HV=NaN(T,nm,nk,nc,nh);HO=HV;HG=HV;
dtv=NaN(T,1); %tracking the time interval in years
for t=2:T
    %Construct daily percentage forward return dF/F 
    %and percentage implied vol change dI/I
    
    %to construct moment estimators, we construct the daily changes backwards over the
    %previous business day
    dt=(ddate(t)-ddate(t-1))/365; %Maturity shrink in years, actual/365 convention
    dtv(t)=dt;
    matdt=mat-dt;
    for c=1:nc
        
        %Today's log(F/S) spread of previous day forward contract,
        %linear interpolation in log spread space
        LFSd=interp1([0;mat], [0;LFS(t,:,c)'],matdt,'linear'); 
        dLFS=LFSd -LFS(t-1,:,c)'; %Daily change in log F/S on previous day contract
        dlf=log(spot(t,c)./spot(t-1,c))+ dLFS; %ln(F(t)/F(t-1))
        df=exp(dlf)-1; %dF/F, percentage change
        
         %First interpolate implied vol across maturity, linear interpolation on total variance
        ivk1=NaN(nm,nk);lfk1=ivk1; lfkp=ivk1;
        for k=1:nk
            dFm(t,:,k,c)=df; %dF/F over previous day, used for moment construction
            
            %Construct the strike of previous day contract
            xo=IV(t-1,:,k,c);tvo=mat.*xo(:).^2;
            lfko=d1(k)*sqrt(tvo)-0.5*tvo; %ln F(t-1)/K
            lfkp(:,k)=lfko+dlf;%lnF(t)/K= ln F(t-1)/K + ln F(t)/F(t-1)
            
            %new vol at d1
            xn=IV(t,:,k,c);  tvn=mat.*xn(:).^2;
            iim=isfinite(tvn); ni=sum(iim);
            tv1d=NaN(nm,1); %only generate numbers on available maturity
            tv1d(iim)=interp1([0;mat(iim)], [0;tvn(iim)],matdt(iim),'linear'); %new vol at d1 at shorter maturity-1d
            ivk1(:,k)=sqrt(tv1d./m1d); %interpolated vol at d1 at mat-1
            lfk1(:,k)=d1(k)*sqrt(tv1d)-0.5*tv1d; %relative strike at d1, LnF(t)/K at d1
            
        end
        %Then interpolate on d1 to get vol at K, spline in Nd1 space
        for m=1:nm
            ym=ivk1(m,:)'; iik=isfinite(ym); ni=sum(iik); %number of available strikes today
            if ni>=3 %at least 3 strikes
                %  IVk(t,m,:,c)=interp1flat(lfk1(m,iik)',ym(iik),lfkp(m,iik)','pchip');
                
                %The following interpolaton on Nd1 should be more stable than on logK because Nd1 is within [0,1]
                ym=ym(iik);
                tvm=m1d(m)*mean(ym).^2; %tvm is juts used for transformation as a scalar
                Nd11=normcdf(((lfk1(m,iik)'+.5*tvm)./sqrt(tvm)));
                
                %as long as we have 3 strikes today, we can extrapolate to get quotes for all
                %five strike from yesterday
                Nd1p=normcdf(((lfkp(m,:)'+.5*tvm)./sqrt(tvm)));
                iiz=isfinite(Nd1p);
                IVk(t,m,iiz,c)=interp1(Nd11,ym,Nd1p(iiz),'pchip','extrap'); %shape preserving cubic spline on Nd1
                %Interpolated implied vol on yesterday's contract
            end
        end
    end
    dIVk(t,:,:,:)=(IVk(t,:,:,:)./IV(t-1,:,:,:))-1; %dI/I over previous day
    
    %Moments on daily changes over fixed number of historical business days
    if t>252
        tt=t-251:t; %information up to time t, changes are all backwards
        x1=dIVk(tt,:,:,:); %dI/I
        x2=dFm(tt,:,:,:); %dF/F
        for h=1:3
            tt=t-hh(h)+1:t;
            x1=dIVk(tt,:,:,:);
            x2=dFm(tt,:,:,:);
            
            nxh=sum(isfinite(x1+x2));
            ii=nxh<hh(h)/2;%more than half data available
            HOt=252*nanmean(x1.^2);  HOt(ii)=NaN;HO(t,:,:,:,h)=HOt;
            HVt=252*nanmean(x2.^2);  HVt(ii)=NaN;HV(t,:,:,:,h)=HVt;
            HGt=252*nanmean(x1.*x2); HGt(ii)=NaN;HG(t,:,:,:,h)=HGt;
        end
    end
    
end

%weekly sampling; every wednesday-- for pricing analysis 
%If holiday, using previous business day information
dv=[ddate(1):ddate(end)]'; 
wdate=dv(weekday(dv)==4); Tw=length(wdate);
dvd=repmat(wdate,1,T)-repmat(ddate',Tw,1);dvd(dvd<0)=NaN;
mdv=min(dvd,[],2); %check gap between the Wednesday and the previous business day if wednesday is a holiday
%double check to see if data gap is large
gap=find((mdv>1))
if length(gap)>0
    datestr(wdate(gap))
end

HVw=interp1(ddate,HV,wdate,'previous'); %return variance
HOw=interp1(ddate,HO,wdate,'previous'); %vol variance
HGw=interp1(ddate,HG,wdate,'previous'); %covariance

Iw=interp1(ddate,IV,wdate,'previous');
Sw=interp1(ddate,spot,wdate,'previous');
LFSw=interp1(ddate,LFS,wdate,'previous');
d1w=interp1(ddate,d1m,wdate,'previous'); %these are just to keep the dimension of the weekly data, no interpolation involved
matw=interp1(ddate,matm,wdate,'previous');


wdaten=busdate(wdate); %The next business day after the pricing date, for initiation of investment
Iwn=interp1(ddate,IV,wdaten,'previous');
Swn=interp1(ddate,spot,wdaten,'previous');
LFSwn=interp1(ddate,LFS,wdaten,'previous');

nobs=sum(sum(sum(sum(isfinite(HVw),2),3),4),5);
figure(1);plot(wdate, nobs,'LineWidth',3);datetick('x','yy')
%starting date, when we start to have moment estimates
T0=find(nobs>0,1,'first')
datestr(wdate(T0))

save('./data/Momentsw.mat','wdate','HVw','HOw','HGw','Iw','LFSw','Sw', ...
    'Iwn','LFSwn','Swn','wdaten','mat','pdelta','currency','nc','nm','nk','cdelta','d1','mdv','Tw','hh','nh','T0', ...
    'd1w','matw','-mat');



